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CN Abstract: PSpectRe is a C++ program that uses Fourier-space pseudo-spectral meth- 

I— I ods to evolve interacting scalar fields in an expanding universe. PSpectRe is optimized 

\^ for the analysis of parametric resonance in the post-inflationary universe and provides an 

w alternative to finite differencing codes, such as Defrost and LatticeEasy. PSpec- 

i-G tRe has both second- (Velocity- Verlet) and fourth-order (Runge-Kutta) time integrators. 

I Given the same number of spatial points and/or momentum modes, PSpectRe is not sig- 

^ nificantly slower than finite differencing codes, despite the need for multiple Fourier trans- 

c/5 forms at each timestep, and exhibits excellent energy conservation. Further, by computing 

I the post-resonance equation of state, we show that in some circumstances PSpectRe 

^_ obtains reliable results while using substantially fewer points than a finite differencing 

^ code. PSpectRe is designed to be easily extended to other problems in early-universe 

^ cosmology, including the generation of gravitational waves during phase transitions and 

Q\ pre- inflationary bubble collisions. Specific applications of this code will be described in 



future work. 
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1. Introduction 

A key challenge for any inflationary model is to ensure that inflation ends and that the 
post-inflationary universe thermalizes, setting the stage for a conventional, hot-big-bang 
cosmology. A promising approach to this problem is provided by preheating and parametric 
resonance, in which the post-inflationary universe undergoes explosive particle production 
[1, 2, 3]. Preheating leads to a fascinating range of phenomenological consequences, possibly 
including superheavy dark matter [4], gravitational- wave production [5, 6, 7, 8, 9, 10, 11] 
and non-Gaussianity [12, 13, 14]. Resonance can occur in a broad swathe of cosmological 
models and, since it depends on the coupling between the inflaton sector and the "rest of 
physics," any observable consequences of resonance constitute a unique window into the 
inflationary epoch. 

Given that resonance is an intrinsically nonlinear process, understanding it is a chal- 
lenging task, and frequently requires numerical studies. Consequently, expanding on a 
theoretical framework initially developed by Khlebnikov and Tkachev [15], a number of 



codes have been developed to study scalar field dynamics in an expanding universe, in- 
cluding LatticeEasy [16] (not to mention the parallelized ClusterEasy [17], and GPU 
accelerated CUDAEASY [18]) and Defrost[19]. These programs all use finite differ- 
encing methods to compute spatial derivatives of the scalar fields. In this paper we de- 
scribe PSpectRe, a Fourier-space spectral code for solving the Klein-Gordon equations of 
two scalar fields in an expanding universe with an arbitrary quartic interaction potential. 
Technically, this is a pseudo-spectral code, in that nonlinear terms in the potential and its 
derivatives are computed by first converting the fields into their position-space representa- 
tion, performing the necessary multiplications, and then taking the inverse transform [20]. 
Since spatial derivatives are taken by multiplying the corresponding Fourier component 
by the wavenumber ki, pseudo-spectral methods are free of differencing noise^, and have 
excellent spatial resolution [20]. On the other hand, overhead from the Fourier transforms 
and inverse transforms exacts a cost in terms of computational efficiency. As currently 
implemented, PSpectRe allows both second- and fourth-order-accurate time integration 
steps. The latter is more computationally intensive, but delivers excellent energy conser- 
vation. PSpectRe is currently implemented to run on a single, shared memory machine, 
but threads efficiently across 16 (effective) cores. The rate limiting step is the fast Fourier 
transform (FFT) algorithm. It would be straightforward to port PSpectRe to a cluster 
environment, given that the FFT libraries we use exist in a fully parallelized form. 

The goal of this paper is to introduce and describe PSpectRe, and to demonstrate 
that it provides a useful cross check on existing finite-differencing-based codes. Looking to 
the future, spectral codes have excellent accuracy and precision, and the investigation of 
preheating-generated observables such as gravitational waves and non-Gaussianity require 
accurate solutions of the background dynamics. The initial version of PSpectRe is equiv- 
alent to LatticeEasy and Defrost in that it solves the Klein-Gordon equation in a rigid, 
expanding background. However, Easther, Giblin and Lim have previously used a spectral 
algorithm to compute the gravitational- wave background generated during resonance [7, 8]. 
This approach naturally lends itself to implementation within PSpectRe, and will be pur- 
sued in future work. Finally, in earlier work Easther and Parry used a pseudo-spectral code 
to study gravitational effects on resonance in a 1+1 dimensional scenario [21, 22]. In terms 
of performance, a "stock" version of Defrost is a few times faster than PSpectRe on 
the same hardware (with a 256'^ grid), but the difference is not dramatic. However, we will 
compute the post-inflationary equation of state, and show that PSpectRe can achieve 
reliable results with a 64^^ grid in circumstances where Defrost requires 128'^ or even 256^ 
points. We attribute this improvement to the lack of differencing noise in spectral codes. 

The format of this paper is as follows. We first review the scalar field evolution equa- 
tions, and the specific program variables used within PSpectRe, which are the same as 
those employed by LatticeEasy. We then describe PSpectRe's basic computational 



^By differencing noise we mean the error introduced wiien computing approximations to a derivative by 
subtracting neighboring-point values and then dividing by some multiple of the grid spacing. Since neigh- 
boring points often have similar values, subtracting them using finite-precision, floating-point arithmetic 
often yields a result with few significant digits. When the result is divided by some multiple of the grid 
spacing, the rounding error is greatly magnified, resulting in a significant source of numerical error. 



engine, and present results on energy conservation and the spatial distribution of energy 
following resonance for the well-studied scenario of an inflation, cp, with a quadratic poten- 
tial, coupled to a massless field, x, via a 0^x^ interaction term. 

2. Equations of Motion 

We consider a homogeneous and isotropic universe with a flat Friedmann- Robertson- Walker 
metric 

ds^ = -df + a^{t){dx^ + dy^ + dz^) . (2.1) 

The scale factor a{t) obeys 

^^ = f., (2.2) 

4vra , , , , 

a =- — (p + 3p). (2.3) 

The universe contains an inflaton field (j){t^ x, y, z) and a generic matter field x(^) x, y, z) 
which obey the Klein-Gordon equation, 

dV •• ■ V^ dV 

U4,+ -=t + ZH4.-^4,+ - = <>. (2.4) 

dV V^ dV 

ox a^ ox 

PSpectRe supports the generic quartic potential 

V = \m^e + \\^4>^ + \m^x^ + Ja^x' + ^s'^'x' • (2.6) 

For this system of interacting scalar fields, the density and pressure are 

P = \{^^ + X') + ^(IV'AP + iVxl') + V , (2.7) 

P = li'P' + X') - ^(IWP + iVxl') - V . (2.8) 

Consider scalar fields whose values and time derivatives are known at the points of an 
N X N X N rectangular lattice, which represents the spatial part of an FRW universe with 
periodic spatial boundary conditions. The same information is contained in the discrete 
spatial Fourier transform, 

c^(fc) = X^0(f)e-^^-, (2.9) 

r 

and similarly for x and X{k). Here, k labels grid points as wave vectors and r labels the 
grid points as spatial positions. The inverse transformation is 



l^^^kyk-r, (2.10) 



k 

The $(/c) are, in general, complex-valued but, since (p and x a-re real valued, (p{r) = (p{r) 
so ^{k) = <&(— fc) and the number of free parameters matches in both representations. 



3. PSpectRe: Algorithm and Implementation 

Working in momentum space, the field evolution is described by A^^ second-order, ordinary 
differential equations for the time dependence of each of ^{k) and X{k), two further degrees 
of freedom for the scale factor (equation 2.3), subject to the overall Hamilton constraint, 
(2.2). Spatial derivatives do not appear, since di4> applied to equation (2.9) simply multi- 
plies $(/c) by fc*, the usual result for the Fourier transform of a derivative. Consequently, 
while finite difference codes work by discretizing the spatial derivatives, (psuedo)-spectral 
methods simply replace V^ with k ■ k. These equations are nonlinear, however, via the 
interaction terms in the potential, and the "friction" terms in the Klein-Gordon equation. 

To avoid confusion, we point out there is a second class of spectral methods for dif- 
ferential equations derived by discretizing the continuum equations and then expanding 
(j){x + 5x) as a discrete Fourier series. This method is used, for example, to solve a linear 
differential equation by transforming the problem into a matrix equation for the Fourier- 
series coefficients [20]. While such methods can produce precision-limited solutions to 
the discretized equations, those solutions are, by definition, only approximate solutions 
to the continuum equations, since the field derivatives have been truncated to some finite 
order. The method used in PSpectRe is distinct from those which use the discretize- 
first methodology, in that we make the transition to Fourier space before discretization. 
This has the benefit of representing the derivative operators exactly, but requires the fol- 
lowing negligibility assumption. Namely, there is a multivalued relationship between the 
discrete-Fourier-series coefficients and the continuum momentum values: a discrete mode 
k corresponds not only to the continuum mode k, but also to the continuum mode k — j^, 
where A^ is the total number of points per grid side and L is the physical length of each 
side of the simulated box. This is because contributions to those two modes are indistin- 
guishable on a finite grid of regularly-spaced samples. PSpectRe uses the convention that 
the first y + 1 Fourier-space components in any dimension represent the modes 0, . . . , ^ 
and the remaining y — 1 points represent the modes ——^ — - through — -£-. By simulat- 
ing systems for which the modes —^ — - through ^^^ — - are negligible, compared to the 
modes — -^^ — through —3-, PSpectRe's representation achieves high fidelity by exactly 
computing the solution for those that are modes explicitly represented. Importantly, the 
dispersion relation induced by the Laplacian operator is not distorted from its continuum 
behavior, as it would be if a discretize-first methodology were employed. 

In the initial phases, resonance typically amplifies a small subset of Fourier modes, 
and in this scenario a Fourier basis is particularly natural. Conversely, one might worry 
that a pseudo-spectral code would be more likely to fail when applied to scenarios in 
which the field has substantial, localized peaks in position space which is possible at late 
times, or in scenarios with topological defects. In work currently in progress, we have used 
PSpectRe to analyze oscillon formation and evolution in three dimensions (c.f. [23]. As 
will be described in a forthcoming paper, our experience is that the pseudo-spectral code 
copes well with very sharply peaked field configurations, even at effective spatial resolutions 
where finite differencing schemes would fail. 



3.1 Field Rescalings and the Velocity Terms 

The time evolution can be handled by any suitable ordinary-differential-equation solver, 
and PSpectRe implements both a second-order (in time) Velocity- Verlet routine [24], and 
a fourth-order Runge-Kutta integrator [25]. To improve the stability of the integration, 
we remove the H(j) and Hx terms from the equations of motion, making use of the same 
rescalings as LatticeEasy.^ We quote them here for convenience 

(j)pr = Aa!'(j) Xpr = -Ao^X ^pr = Bx dtpr = Ba^dt . (3-1) 

If A^ is 0, A;^, is 0, and m^ is nonzero, then 

B = m^ r = 3/2 s = 0. (3.2) 

If A^ 7^ and A;^ is 0, we set 

B = y%4>o r = 1 s = -l, (3.3) 

where (pQ is the initial value of the amplitude of the k = mode for (p. In all cases, 

A = ^. (3.4) 

since the (j) field is initially dominant. We also follow the LatticeEasy convention of 
using a prime to denote differentiation with respect to the rescaled program time, while 
dots denote differentiation with respect to physical time. Rescaled (or "program" ) variables 
are written with a subscript "pr" . 

3.2 Transformed Equations of Motion 

To write the equations of motion for the Fourier coefficients of the fields, we take the 
Fourier transform of the Klein-Gordon equation. This results in the following equation for 
the amplitude of each momentum mode with wave vector k : 

$ + ^c,_3^c, + § = o, (3.5) 

where the hat above the partial derivative indicates a Fourier transform. The equation for 
X{k) is of precisely the same form. 

In terms of rescaled variables, the equations of motion are 

$;; = -{s - 2r + 3)-$;, + Hs - r - 2) - + —]%r 

- \kpr\'a-''-^<^pr - ^ , (3.6) 



where 



dVpr a 



bpr 52 



upr 



2s— 2r \ 2 



'V' ^ j^2 pr ^ j^2 pr^Pr 



(3.7) 



■^A full explanation of these rescalings can be found in the LatticeEasy documentation. 



The expression for -qj^ is similar, with $pj. and Xp^ exchanging roles. The variable 
rescalings are designed so that s — 2r + 3 = 0, so the <l>p^ dependence is eliminated. 

The scale factor a is evolved using the same version of the Friedmann equations em- 
ployed by LatticeEasy: 

a = -2^ + 87ra(^ (|V0|2 + |Vxp) + aV) (3.8) 

where angle brackets denote averages over the simulation box. In terms of rescaled vari- 
ables, this becomes 

a" = {-s - 2)^ + ^a-^'-^'-'(^{\VpApr? + WprXpr?) + a'^+'v) . (3.9) 

3.3 Nonlinear Terms 

Spectral codes are simple to implement in linear systems, because the individual momen- 
tum modes do not mix. However, preheating and resonance are driven by the nonlinear 
couplings between modes, and properly accounting for these is the primary challenge faced 
by any spectral treatment of resonance, both in terms of algorithmic complexity and exe- 
cution time. 

We have to compute several different nonlinear terms. Firstly, to evolve the scale 
factor we need averages over the entire box for both quadratic (e.g. (|V(;^p)) expressions 
and the potential, which - in principle - can contain mixtures of arbitrary order in the 
fields. Writing a variable (e.g. (j)) in terms of its Fourier expansion and integrating over 
the box generates delta functions, and we evaluate the corresponding terms as a sum of 
squares: 



box 



fc_spacc 

/ l^'^l' = 4 E l^l'l'^(^)l'- (3-10) 

./box -' ' ^ 

fc_space 

Integrating the potential is trickier, since we now have cubic and higher-order terms. In 
principle, we could evaluate an M — th order term by writing ^{ki)^{k2) ■ ■ ■ <I>(/cjvf); per- 
forming the spatial integral gives a delta function which ensures that fci + /c2 + • • ■ + kM = 0, 
and we could then sum over all the corresponding momenta that satisfied this expression. 
However, the number of possible combinations is prohibitively expensive as M becomes 
large. Consequently, we shift the fields into the position space representation, compute 
the potential, and integrate over the (spatial) box using Simpson's rule. This approach 
is necessarily inaccurate, since the product of the $(/ci)<I>(A;3) • • • $(fcM) will contain terms 
of frequency up to ^ Mkmax, which will not be resolved on our spatial grid. In other 
words, if <I> contains a Fourier-space component e """-^^j then <I> will contain a contribu- 
tion (^e^f'maxx'j _ Qi{Mkmax)x^ jj^ practlcc, power in such off-grid modes is automatically 
aliased by the DFT to some mode represented on the grid. PSpectRe contains an option 
to "pad" the spatial grid for the computation of the energy integral, so that it contains 



{pN)^ points, where p is a (small) integer. The padding and unpadding algorithms are 
specified in the appendix. This slows the code, but dramatically improves instantaneous 
energy conservation, as we discuss below. 

To check the energy conservation we need to compute the average density, (p) over 
the simulation volume. Quadratic terms are simply dealt with via Parseval's Theorem. In 
terms of rescaled variables, p is given by 



^{(Ppr"^ + Xpr^) - r — {(ppr(l)pr + XprX'pr) + 2'^^ \~ 



+ Q, y\\' prYprl ~r \^ prXpr\ ) 

+ V)- (3-11) 

We use four applications per field of Parseval's theorem to compute every term in the 
energy density except for the potential: 



fc-space 
fe-space 

'pr4'pr = ^ Yl [Re( V)Re($;,) + Im( V)Im(^;,)] 

box -' ' ^ 

fc-spacc 
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I v pr(Ppr\ 



-^ Y, IVI'IVP. (3.12) 

fe-space 



Finally, we need the Fourier components of the potential derivatives, ^ and ^, which 

are generally nonlinear with respect to the fields, including terms proportional to x'^4': 0^X; 
4>^ and x^- We handle these terms by transforming the fields into position space, performing 
the multiplication (such as x^*/*)) s-^id then transforming the result back to Fourier space to 
proceed with the spectral evolution. This is necessarily time-consuming, and is the price 
we pay for evolving a nonlinear system of equations with a spectral algorithm. 

3.4 Integrators 

PSpectRe provides two user-selectable integrators: A standard fourth-order- in-time Runge- 
Kutta integrator and a second-order-in-time Velocity- Verlet integrator. The Runge-Kutta 
scheme is more expensive, since the nonlinear terms in the potential must be transformed 
four times per time step. The Velocity- Verlet scheme is less expensive, but only second- 
order accurate in time. The scale factor is evolved along with the fields using the same 
integration scheme. In the case of the Velocity- Verlet scheme, this violates one of the 
scheme's assumptions, that the second derivative does not directly depend on the first 
derivative, since a explicitly depends on d. Since a has become small and slowly-varying 



by the time resonance sets in, the effect of this discrepancy is minor, but reduces the Verlet 
integrator to first order accuracy for the scalar factor. 

3.5 Initial conditions 

PSpectRe provides two routines for initializing the fields, which emulate the procedures 
of LatticeEasy and Defrost, respectively. In either case, the initial value for the 
k = mode is provided by the user in physical units, and sets to the mean value of the 
corresponding field. In our experience, the post-resonance state of the systems we have 
simulated are qualitatively, if not quantitatively, unaffected by the choice of the initial- 
conditions routine. 

3.5.1 LatticeEasy Emulation 

We refer the reader to the LatticeEasy documentation (section 6.3) for a complete dis- 
cussion of the assumptions and approximations that underlie the following formulae. The 
initial values of all mode amplitudes except the k = mode are assumed to be given by 
Rayleigh distributions with RMS values W^ such that 



where 



N^)]' = l^l' + S- (3.14) 



In terms of rescaled variables, the RMS values are 

ly.^^ = ^ , (3.15) 

where 

[ojprikpr)] = \kpr\'^ + J"' . (3.16) 

Once we have computed these RMS amphtudes, we generate a uniform deviate Y and use 
it to compute a randomized quantity fr : 

Finally, we generate two random phases, 9i and Or, for the left-moving and right-moving 
components of the momentum mode, and set the real and imaginary parts of its amplitude 
via 

These initial values for the momentum-mode amplitudes must be multiplied by a factor of 
N'^ in order to be consistent with the Fourier transform conventions listed above. Since 
the fields are real-valued, we must ensure that these initial momentum mode amplitudes 
obey the appropriate conjugate symmetry. 



For simplicity, the initial values of the time-derivatives of the k = mode amplitudes 
are chosen to be zero for both fields. The initial values of the derivatives for the other 
modes are set according to the formula 



$'^ = ^j;[± ioJpr{kpr) + {r- l)Hpr\ , 



(3.19) 



where Hpr = ^. 

3.5.2 Defrost Emulation 

Defrost uses a convolution-based method to initialize the fields and their derivatives. 
The initialization process is viewed as constructing a sample from a Gaussian random field 
defined as 



ip{x,t) 



d^k 



Akx 



bk cos cot + Cfc sin cot 



(27r)5 J V2uj 
where the mode amplitudes are uncorrelated random variables: 



and 



The effective mass is 



ihK') = {ckcl') = 5{k- k') 
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(3.20) 

(3.21) 

(3.22) 
(3.23) 



This is implemented by convolving white noise with a (spherically symmetric) kernel 
function 



i{r) 



1 



(2^)1 



Sk 



ikx 



2uj 
k'^dk 



(3.24) 



sin fcr 



(A;2 + mlo) 4 ^'f' 



Unfortunately, this kernel function is singular as r — )■ 0, so it is regularized by introducing 
a Gaussian cut-off at some scale q below the Nyquist frequency 



e(r) 



1 



k dk sin kr 
(A;2 + m^fj) 4 kr 



exp 



n2 



(3.25) 



Unhke Defrost, which uses a cut-off of ^^Nyquist, PSpectRe uses a cut-off of ^/^Nyquist- 
The larger cut-off produces field-mode amplitudes closer to their physical values. The 
regularized kernel function, ^(r), is then convolved with the {bk}: 



^{x, 0) 



1 



(27r)S 



d''kd-'x'bk^{x')e 



r^\Jk{x-x') 



(3.26) 



A similar procedure is used to initialize the velocity fields. 
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Figure 1: We plot the energy conservation for two 32"^ runs for the quadratic model. The red 
line uses the Verlet integrator, blue shows the Runge-Kutta results. Before resonance (the region 
shown in the lower plot), the Runge-Kutta integrator is far more accurate, but both models see 
significant departures from exact conservation once resonance begins. For both runs L = 2 and the 
timestep is 0.005. 

3.6 Units 

The masses and coupling-constants appearing in the potential and the field amplitudes 
are specified in physical units normalized such that the Planck mass is 1. The size of 
the box and the program time are specified in terms of the (rescaled) program units. For 
a potential dominated by an inflaton mass term, this is equivalent to specifying these 
quantities in terms of the inflaton mass. Specifically: 

L = ^ . (3.27) 

4. Validation 

PSpectRe was validated by comparing to LatticeEasy and Defrost using the m?(jp'- 
model, and also by ensuring that its own results are self-consistent as the timestep, inte- 



10 



0.00015 : 
0.00010 : 
0.00005 - 



"yvv^vv<;vvOvv^vvOwvV/wl/vv\>vv^wOw 



50 



100 



-0.00005 
-0.00010 
-0.00015 




Figure 2: We plot the same scenario shown m Figure 1, with a padding factor of 2. The energy 
conservation improves by an order of magnitude. 



grator and box size are varied. We primarily focus upon the quadratic model 



(4.1) 



which is the default model in Defrost, and all runs shown here have the same parameter 
values as those presented in [19]. We also tested our code against the quartic A(/>^ model 



^(<A,X) = ^5V/ + ^A<^^ 



(4.2) 



which is implemented inside LatticeEasy. Formally, there is no conserved energy in the 
system, but we gauge the accuracy of the evolution by checking the accuracy with which 
the averaged Friedmann equation is satisfied, plotting 



(P) 



5p 



1, 



(4.3) 
, initially ~ 10~^) when energy is "perfectly" 



a quantity that will vanish (up to 
conserved. 

Figure 1 shows the energy as a function of the integration scheme. Initially the Runge- 
Kutta greatly outperforms the Verlet, but both degrade as resonance hits. We compute 
the average potential in our box by constructing y((/>, x) is real space and performing a 
numerical integral. However, for this case the shortest momentum modes have a non- 
trivial amplitude and the resulting integral effectively under-samples the "fine structure" 
in the potential. To test this hypothesis we have added a "padding" feature to the code, 
where the A^^ momentum modes are Fourier transformed back into a (pN)^ spatial lattice, 
which effectively interpolates extra points into the numerical integral. The padding and 
unpadding algorithms are specified in the appendix. Figure 2 shows the same scenarios as 
those plotted in Figure 1 with a padding factor p = 2. In this case, the energy conservation 
for both integrators is improved by an order of magnitude. Finally, Figure 3 shows 64^ 
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Figure 3: We plot the same scenario shown in Figure 1, for Runge-Kutta integrators and a box 
size of 64. The red plot is unpadded, and the blue plot has a padding factor of 2: we show the 
latter simulation on its own in the lower panel. In this case there is an initial transient which does 
not depend on the algorithm, and the subsequent evolution is followed at better than 1 part in 10^. 



runs with and without padding. With padding we obtain an accuracy of the order of parts 
in 10^ throughout the run. Increasing A^ (for constant box size, L) will often have a similar 
effect to using padded grid for the potential, if the higher momentum modes have a low 
amplitude throughout the simulation. However, using the padded grid increases the run 
time by a factor of approximately 2.5, whereas doubling N will increase runtime by an 
order of magnitude, since the higher resolution simulation tracks eight (= 2'^) times as 
many momentum modes as the lower resolution simulation. 

More generally, our results here reflect several paradoxes associated with the discussion 
of energy conservation in a preheating code. Firstly, these systems are intrinsically non- 
linear, and mode-mode couplings transfer power to momenta that are unresolved within 
our simulation. Consequently, a code that has perfect energy conservation while tracking 
a finite number of momentum modes is unphysical to the extent that it suppresses this 
transfer of power, and effectively imposes a low-pass filter on the dynamics. Our use of 
a padded grid avoids this problem, since we can accurately compute the potential-energy 
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Figure 4: Energy conservation performance of Defrost [blue] and PSpectRe [red]. Both codes 
are run with 256'^ points and L — 10 for Defrost's defauh quadratic model. 



integral. Note that we evaluate the quadratic expressions for the momentum and gradients 
exactly without resorting to padding - writing these spatial integrals with the appropriate 
Fourier expansions creates a sequence of S functions and our integrals are evaluated by 
summing the squared-amplitudes of the Fourier coefficients. Finally, this code (like most 
other preheating simulations, see [28, 29] for exceptions) implicitly ignores any contribu- 
tion from metric perturbations to the scalar field evolution. Despite these caveats, we have 
demonstrated that PSpectRe can conserve the Hamiltonian constraint for this system to 
high accuracy, if desired. In practice, PSpectRe's Verlet integrator is suitable in most 
situations, and offers savings in both runtime and memory use. On the other hand, if high 
accuracy is needed, the Runge-Kutta integrator provides better performance than setting 
the Verlet timestep to a very small value. 

We compare the energy conservation of PSpectRe and Defrost in Figure 4. We find 
that we can make the timestep in PSpectRe somewhat larger than Defrost's "default" 
value and obtain a comparable level of energy conservation. For this case Defrost ran 
approximately twice as fast as PSpectRe, with both codes saturating a 16 (effective) 
core, shared-memory machine. In practice, performance can degrade substantially when 
"additional" variables are computed, such as the Newtonian gravitational potential, or 
the user calls for frequent data output. Given that the Fast Fourier Transform is an 
0{3N^logN) process we would expect PSpectRe to run more slowly than a comparable 
finite differencing code, which naively scales as 0{N^), but in practice the time penalty is 
not exorbitant. 

While PSpectRe evolves the fields in momentum space, it can output its results 
in terms of the spatial values of fields, and of the various classes of energy within the 
model. As an example. Figure 5 shows the output of a quadratic potential with L = 3 
at 256^, and we recover the spatial patterns described in [27]. Parametric resonance is 
typically associated with a fixed physical scale, which determines the length scales that 
the simulation must accurately resolve. Clearly, algorithms which allow this characteristic 
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Figure 5: The top plot shows the relative energy density on a two dimensional slice through the 
box at the onset of resonance in a quadratic model; the lower plot shows contours at which the 
overall density is 1.5, in units where the average density is unity. These plots are made with L = 3 
and 256'^ spatial points. 

scale to be close to the smallest length (or highest momentum) are more memory-efficient 
than those which require a substantial amount of separation between these two scales. As 
noted in the introduction, spectral codes compute spatial derivatives using an algebraic 
identity, rather than by comparing values of the fields at adjacent points, and are thus 
not affected by differencing errors. Consequently, it is reasonable to hope that PSpectRe 
might need fewer points (or momentum modes) than a finite-differencing code applied to 
the same system. To investigate this, we extract the instantaneous equation of state w as 
a function of time [19], where 



w 



M 
(p) 



(4.4) 



In Figure 6 we show the computed value of w for the quadratic model with box size L = 
5, using Defrost at A^ = 128 and N = 256. The two curves almost completely overlap, 
showing that the code is effectively in the continuum limit when A > 128. In Figure 7 we 
plot the same quantify at smaller N, and find that at A = 64 and A = 32 the resulting 
values for w differ substantially. Conversely, with PSpectRe we find excellent agreement 
between A = 128 and A = 64, while the A^ = 32 case is still largely reliable. Consequently, 
it appears that pseudo-spectral code does allow us to accurately follow the dynamics of 
resonance with a lower cutoff length scale. From a practical perspective, this can either 
have a substantial effect on the runtime, given that at a constant timestep, performance 
naively scales as A^^, or (alternatively) PSpectRe can simulate larger spatial volumes 
than a finite differencing-code, given a fixed amount of computational resources. Finally, 
while PSpectRe and Defrost both reach a "continuum limit" when the number of points 
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Figure 6: We show w, as a function of time, in a Defrost simulation at L = 5, with all other 
variables at their default values. The black line is a 256'^ simulation, and red points are from an 
otherwise identical computation at 128^. 
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Figure 7: The left panel shows w three Defrost simulations, at 32"^ (red), 64^ (blue) and 256^ 
(black, solid), whereas the right panel shows the results of three PSpectRe simulations, at 32^ 
(red), 64'^ (blue) and 128'^ (black). Other than the number of points, the parameters are identical 
between these runs and those in Figure 6. Note the vertical scales are different between the two 
plots, and the time interval shown is subset of that used in Figure 6. 

being tracked in the simulation is sufficiently large, the two codes do not overlap with each 
other as the detailed time-evolution of w depends very sensitively on the implementation 
of the initial conditions. However, we do recover the late time limiting value w ~ 0.2 seen 
by Frolov. 



15 



5. Conclusion 

We have demonstrated that PSpectRe yields a highly-accurate and computationally- 
efficient approach to the solution of scalar-field evolution during parametric resonance. 
PSpectRe thus provides a useful contrast to finite-differencing algorithms such as those 
in Defrost and LatticeEasy. 

The present discussion is intended primarily as a "proof of concept" , demonstrating 
that a pseudo-spectral algorithm is a practical method for the analysis of preheating, de- 
spite the extra computational overhead introduced by the nonlinear terms in the evolution 
equations. In doing so, we have recovered a variety of existing results in the literature. 
We find that PSpectRe's energy conservation performance is excellent, and scales in well- 
understood ways as the grid size and timestep are varied. Moreover, we demonstrate that 
some energy is "lost" to the simulation when it is transferred to high frequency modes. 
While PSpectRe does not evolve these modes, it can account for them by "padding" the 
spatial integral used to compute the potential energy. 

As a pseudo-spectral code, PSpectRe does not rely on comparing values of fields at 
adjacent lattice points to take spatial derivatives. Consequently, in some circumstances 
PSpectRe provides accurate results with smaller grids than those needed by finite differ- 
encing codes. When the grid sizes are equal, PSpectRe is still competitive with finite- 
differencing codes in terms of runtime, but the repeated Fourier transforms do introduce a 
noticeable performance penalty. Reducing the grid size by a factor of 2 naively improves 
runtime by (slightly more than) a factor of 8, with fixed timestep, while maintaining the 
same effective resolution as a finite-differencing code. Consequently, by changing runtime 
settings, PSpectRe can be optimized for precision, speed, or the ability to run simulations 
over large spatial volumes. 



A. Download and License 

PSpectRe is publicly available^, and is made available under a BSD-style license. Its 
current features are on a par with those of Defrost and LatticeEasy, and we plan to 
extend its capabilities in future work. PSpectRe compiles with both the freely-available 
GNU C-|--|- compiler and the Intel C-|— |- compiler. The code calls either the FFTW library, 
version 3 [26] or Intel's Math Kernel Library to compute the Fourier transforms Our experi- 
ence is that the Intel library outperformed FFTW on machines where both were available. 
If an extended precision FFTW library is available^, the user may select a mode in which 
all computation is done using extended floating-point precision. PSpectRe has been run 
on both Mac OS and Linux based systems. Finally, we warn that not all combinations of 
features and settings have been exhaustively tested. 

See: http : //easther . physics . yale . edu/ 
*The extended precision build is often referred to as the "long double" build after the name of the C 
extended-precision data type. 
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B. Padding and Unpadding Algorithms 

The conjugate symmetry allows the DFT of a grid oi N x N x N real values to be stored 
using a Fourier-space grid of only A^xA^x(-2- + l) complex values. These values are further 
restricted such that they contain exactly A^^ degrees of freedom. The padding algorithm, 
which interpolates between a Fourier-space grid of size N x N x (y + l) and one of size 
pN X pN X pN for some positive integer p, is specified in Algorithm 1. The unpadding 
algorithm, which reverses the procedure, is specified in Algorithm 2. 
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Algorithm 1 Padding Algorithm. In the actual implementation, the input and output 
grids use the same underlying memory buffer. 



Define: S{v) 



v>pN- ^ 
otherwise 



v-{p- 1) A 
undefined 



\ v = 

I otherwise pN — v 
Require: input[0 ... (A - 1)] [0 . . 
Require: output[0 . . . {pN - 1)][0 
for X = pN — 1 to do 
for y = pN — 1 to do 
for z 



(A-i; 

. . . (pA 



][0...f] 
-1)][0.. 



2 J 



N 



then 



output[x][y][z] 



2^ to do 
if S{x) is defined and S{y) is defined and z < ^ 

(z = f '^inpnt[S{x)][S{y)][z] 

1 otherwise input[S'(x)][S'(y)][z] 
else 

output[x][y][z] -^ 
end if 
end for 
end for 
end for 

for a; = to pN do 
for y = to pN do 

z •(— {Set the z = (and z = ^) modes by conjugate symmetry.} 
if S{x) is defined and S{y) is defined then 

if S{xc) is not defined or S{yc) is not defined then 
if (x, y) / (xc, yc) then 



ioutput[a;][y][2;] 
— output[x][y][z] 



output [x] [y] [z] -f 

output [a;c][yc][^] 
end if 
else if S{xc) is defined and S{yc) is defined then 

Assert: output [xc][yc] [-2] = output [x][y][z]* 
end if 
end if 
end for 
end for 
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Algorithm 2 Unpadding Algorithm. In the actual implementation, the input and output 
grids use the same underlying memory buffer. 



Define: S{v) 



v>pN- f v-{p-l)N 



otherwise undefined 

{ v = Q 

I otherwise pN — v 
Require: input[0 . . . {pN - 1)][0 . . . {pN - 1)][0 . . . ^] 
Require: output[0 ... (iV - 1)] [0 ... (A^ - 1)] [0 ... f] 
for a; = to pN do 
for y = to pN do 

z -^ {Restore the z = modes.} 

if <S'(2;) is defined and S{y) is defined then 

if S{xc) is not defined or S{yc) is not defined then 
if {x,y) / {xcVc) then 

input [x] [y] [z] ■(— 2 input [x][y][z] 

input [xc][yc]N ^ 
end if 
end if 
end if 
end for 
end for 

for a; = to pN — 1 do 
for y = to pN — 1 do 
for z = to ^ do 

if S{x) is defined and S{y) is defined and z < ^ then 

output[5Wl[5fe)|H^(^=*. 2"P"'H[!'1W 

I otherwise input [xj [yj [zj 

end if 
end for 
end for 
end for 
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